Linear Stability of Scroll Waves 



Hcrve Henry and Vincent Hakim 
Laboratoire de Physique Statistique, Ecole Normale Superieure, 
24 rue Lhomond, 75231 Paris Cedex 05, France 
(February 8, 2008) 



O 
O 
O 

(N 

> 
O 

in 



in 
d 



> 
o 
m 
o 



o 
o 



> 

x 



A full linear stability of a straight scroll wave in an ex- 
citable medium is presented. The five eigenmode branches 
which correspond to deformation in the third dimension of 
the five main modes of two-dimensional (2D) spiral dynamics 
are found to play a dominant role. Modulations in the third 
dimension have stabilizing or destabilizing effects on the dif- 
ferent modes depending on the parameter regimes. For un- 
twisted scroll waves, our numerical results confirm the rela- 
tion between the long-wavelength behavior of the translation 
branches and 2D spiral drift in an external field but show 
no similar direct relation for the meander branches. The in- 
fluence of twist on the different branches is investigated. In 
particular, the sproing instability is seen to arise from the 
twist induced deformation of the translation branches above 
a threshold twist. 



Scroll waves are three-dimensional (3D) extensions of 
the familiar spirals |t) of excitable media. They have been 
directly observed in three-dimensional chemical reactions 
in gels and are suspected to play a role in the slug 
phase of slime molds life cycle || and most importantly 
in ventricular fibrillation j|]5). Numerous simulations 
of scroll waves have been performed, different instabil- 
ities have been noted and a somewhat similar phe- 
nomenology has been found for 3D complex Ginzburg- 
Landau vortices ||. However, even in the simplest case 
of an homogeneous and isotropic excitable medium, the 
mechanisms of the different instabilities, their system- 
atics and the influence of scroll wave twist still appear 
poorly understood. We attempt here to bring some clar- 
ification by performing a full linear stability analysis of 
straight scroll waves in various parameter regimes. In 
contrast to earlier stability studies [p|jlc|], this gives ac- 
cess not only to the most unstable eigenmode but also to 
the other, stable or unstable, modes. In particular, this 
enables us to follow the five branches of modes which 
emerges from the two translation modes, the rotation 
zero mode and the two Hopf bifurcation meander modes 
of 2D spiral motion. As for spirals |[lj, these are found 
to play the dominant role in scroll wave dynamics. 

The excitable medium dynamics is described in a 
simplified but usual way by the two-variables reaction- 
diffusion system 



d t u = V 2 u + f(u,v)/e 
dtv = <5V 2 v + g(u, v) 



(1) 

(2) 



We restrict ourselves to the singly diffusive case with 5 = 
and choose for definiteness Barkley's reaction terms 
fll2| with f(u, v) — u(l — u)[u — (v + b)/a], g(u, v) = u — 
v which permits fast direct simulations and comparison 
with some previous results. Previous experience with 
2D spirals lead us to expect that similar results would be 
obtained with a different choice of excitable kinetics. The 
results should also be, at least qualitatively, comparable 
to experiments using Belousov-Zhabotinsky reaction [fl3| . 

Our numerical method is analogous to the one used in 
a previous linear stability analysis of spiral waves [p"l|| . 
We introduce a rotating coordinate system twisted along 
the z-axis where the steady scroll solution is time in- 
dependent. That is, we seek u and v as functions of 
r,4> = 9 — u)t — T w z,t and z where (r, 9, z) are cylindri- 
cal coordinates and t w is the imposed twist. In these 
coordinates, Eq. ([I]j2]) read 

(d t + 2t w 8 2 z - d 2 zz )u = (woV + r 2 d% + V 2 2D )u + f(u, v)/e 
d t v = loO^v + g(u, v) (3) 

The Laplacian is two dimensional on the r.h.s of Eq. (||) 
[in (r, <p)] in contrast to the three dimensional one in 
Eq. (0). The twist t w can be chosen at will and is not 
constrained by the finite size of the simulation box in the 
z-direction as in usual direct simulations. Determining a 
steady scroll rotating at the frequency u) = uj\ consists 
in finding («o(r, (/)), vo(r, </>), uji), that nullifies the r.h.s 
of Eq. (||). To solve this nonlinear eigenvalue problem, 
the r.h.s. of Eq. (|3|) is discretized on a Nr x finite, 
2tt 0— periodic box of radius R using finite differences 
of 4 th order in space for the differential operators with 
rotationally symmetric boundary conditions imposed at 
the box edge, d r uo\ r =R — 0. The non-uniqueness of 
uo(r,4>), vo(r, (f>) brought by the rotational invariance of 
(^) is taken care of by arbitrarily setting the value of 
^(Nj-jN^) equal to 0.5. Therefore the discrete problem 
is reduced to finding the zero of a (2 x N r X iV^) valued 
function of the (2 x N r x N^ — l) unknown values of uo and 
vq and u>\. This is solved by Newton's method starting 
from an initial guess of the solution provided by a direct 
simulation of the evolution equations. In a box of size 
60 x 120, solutions of the discrete problem are obtained 
with an accuracy of I0~ 8 (in a L2 norm). 

Once a steady scroll wave (uq, vo 7 u>i) is found, the time 
evolution equations (|^) are linearized around (uq,vq) in 
the twisted rotating frame. Setting u = uo + exp[er(fc z )i — 
ik z z]u\{r, 4>), v = vq + exp[a(k z )t — ik z z]v\{r, 4>), leads to 
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a linear eigenvalue problem for the growth rates a(k z ) as 
a function of the wave- vector k z , 

a(k z ) u x = (-kl + 2iT w k z d 4 ,)u 1 + (ui^ + + X7 2 2D ) 

+ [d u f(u ,v )ui + d v f(u Q ,v Q )vi}/e 
a(k z ) vi = Ldid^vi + [d u g(u Q , v )ui + d v g(u ,v )vi] 

These linear equations are discretizcd in the same man- 
ner as Eq. (||) . We denote by C kz the discrete version of 
the linear operator appearing on the r.h.s. of Eq. (|J). 
Being interested in the stability of the scroll wave and 
in its modes of destabilization, we focus on determining 
the eigenvalues of C k , of largest real parts. This can be 
accurately done by using an iterative method fully de- 
scribed in juj. Firstly, the contributions of eigenvalues 
with very negative real parts are effectively suppressed by 
computing by iterations x\ = (1 + dtCk z ) to ^ dt xo, for an 
arbitrary vector xq and a sufficiently large integer to/dt 
and a sufficiently small dt (in the following calculations, 
typical values are to = 5 and dt ~ 10 ~ 5 ). One then 
builds an orthonormal set of vectors that spans {aq, x<i = 
(1 + dtC k J^ dt x , x m = (1 + diAJ (m_1) * lM zo} 
and evaluates the restriction of Ck z to this sub- 

space. The integer t\jdt should be chosen large enough 
to make (1 + Ck z dt) tl / dt significantly different from the 
identity but small enough to limit computation time (as 
a typical value, t\ = 0.5 is used here). One finally 
computes the n eigenvalues {<7j,l < j < n} of largest 
real part (Re(<7j) > Re(<Jk),j < k) and the correspond- 
ing eigenvectors of . They approximate the corre- 
sponding n eigenvectors of C kz with an accuracy of order 
exp[— (<t^ — cr^ 1 )(m — n)to]. With m = 50, we reach a 
precision of more than 10 -6 , in computing eigenvectors 
of C kz (in La norm \\[C kz - a(k z )](u u vx)\\ < 10" 6 ). 

We have determined steady rotating scroll waves and 
computed their stability spectrum for several different 
parameters of Eq.([j],||) where steady spiral waves exist 
and are stable. We now summarize these results begin- 
ning with the case of untwisted scroll waves. 

For 2D-spiral waves, the linear stability operator has 
three marginally stable modes coming from the symme- 
tries of the dynamics : the rotation eigenmode with an 
eigenvalue equal to zero and the two translation eigen- 
modes with purely imaginary eigenvalues equal to Hut, 
the spiral rotation frequency. For untwisted scroll waves, 
these correspond to eigenmodes of the linearized operator 
at k z — 0. The computations reported in Figs. [j],|]J^ show 
good agreement with these expectations. For k z — 0, we 
find an eigenvalue of magnitude smaller than 10~ 9 which 
we identify with the rotation eigenmode and two eigen- 
values equal to Huii with an error of 10 which we 
identify with the translation eigenmodes |lq ] . 

For 2D-spiral waves, in addition to these three 
marginal eigenvalues, two complex conjugate eigenval- 
ues cross the imaginary axis on the 2D meander insta- 
bility line. The corresponding meander eigenmodes play 
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FIG. 1. (a) Real and (b) imaginary parts of the domi- 
nant eigenvalues of the spectrum as functions of k z : trans- 
lation (+), rotation (•) and meander(o) branches. The ki- 
netics parameters are a — 1., b = .01 and e = .05 and 
qualitatively give 2D spirals with small core far from the 
meander line. No instability is observed. The spiral rota- 
tion frequency is wi = 1.103 and the width of the excited 
zone between two u = .5 iso- lines far from the tip corre- 
sponds to a wavenumber k\ ~ 0.90. A fit for k z <C 1 of 
the translation branch gives at — 1.1 i+ (—0.38 + 1.31 i) kl in 
agreement with the independently measured drift coefficients 
aii = 0.37, a± — 1.32. A similar fit for the meander branch 
gives a m ~ -0.18 + 1.05 i + (-0.49 - 0.49 i) k\ 



an important role in spiral dynamics and belong to the 
untwisted scroll wave spectrum at k z = 0. The results 
of our computations show that in general, the eigenval- 
ues with the largest real parts in the untwisted scroll 
wave spectrum belong to the five finite-fc 2 branches of 
modes originating from the translation, rotation or me- 
ander eigenmodes at k z = 0. The rotation mode appear 
to be stabilized by finite k z values. On the contrary, we 
have observed three different behaviors of the translation 
and meander branches as a function of k z for different 
choices of the parameters a, b, s characterizing the reac- 
tion term f(u,v)/e in Eq. (|l|): 

-In the first case, both the translations and the mean- 
der modes are stabilized when k z becomes non-zero, as 
shown in Fig. [|. 

The second case is shown in Fig. ^. The translation 
modes now become unstable for low k z while the mean- 
der modes are stabilized. This corresponds to the "nega- 
tive filament tension" instability previously described in 

0- 

- A third case opposite to the previous one is also ob- 
served as shown in Fig. ||. The translations modes are 
here stabilized when k z becomes non-zero. On the con- 
trary, the real part of the meander modes grows pro- 
portionally to k z at low k z before decreasing for higher 
wave numbers. This qualitatively corresponds to the phe- 
nomenon reported in the scroll waves can be unsta- 
ble with respect to meander at finite k z in a parameter 
regime where 2D spirals are stable, as shown in Fig. [|b. 

The fourth logical possibility consisting of a destabi- 
lization of both the translation and meander modes was 
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FIG. 2. (a) Real and (b) imaginary parts of the domi- 
nant eigenvalues of the spectrum as functions of k z : trans- 
lation (+) and rotation (•) branches. The kinetics pa- 
rameters are a = .96, b = .001 and e = .1 and qual- 
itatively give a weakly excitable medium with large core 
non-meandering spirals. The meander modes are not plot- 
ted because their real parts is more negative than the values 
shown. The translations branches are unstable for low k z . 
With the notations of Fig. |, a t ~ .52 % + (0.88 + 1.26 i)k\ 
for k z <C 1 and a\\ = — 0.87, a± — 1.25. Other values are 
wi = 0.521, k x ~ 0.95. 



not observed in the present limited search. 

We have performed some direct numerical simulations 
of Eq. (|l],||) to observe the nonlinear development of the 
different instabilities p2[ . When the translation modes 
are destabilized at low k z (Fig. |J), the filament grows 
and no restabilization is observed. On the contrary, a 
restabilization at finite amplitude is observed when the 
meander modes arc destabilized at finite k z . This 3D 
instability thus inherits the direct character of the 2D 
meander bifurcation. These observations corroborates 
previous findings (7|^]. However, a detailed characteri- 
zation of the nonlinear dynamics is beyond the scope of 
the present letter. 

We now proceed and study the effect of a finite 
twist (t w 0) |[(|. At a general level, we note that 
the spectra obey cr(— k z ) = <r*(kz) where the star de- 
notes complex conjugation, since £_ ^ = C*_ k . As 
in the untwisted case, the symmetries of the dynam- 
ics provide three known eigenvalues, one zero eigen- 
value at k z = for the rotation eigenmode and two 
complex conjugate eigenvalues for the translation eigen- 
modes: iuj\ at k z = —t w corresponding to the eigenvec- 
tor exp(i(j))(d r UQ + id^uo/r, d r VQ + id^v^/r) and — iw\ at 
k z = +t w for the complex conjugate eigenvector. These 
exact results provide sensitive checks on the numerics. 

We generally find that increasing the twist r w from zero 
reduces the stability of the translation branches, destabi- 
lize them or amplify the instability when they are already 
unstable. The effect of twist increase appears much less 
pronounced on the meander modes. We focus here on 
the most interesting parameter regimes such as those of 
Fig. [j],||a for which the untwisted scroll waves are stable. 
A representative case is shown on Fig. 0. The kinetics 




0.25 0.5 0.75 1 " "0 0.25 0.5 0.75 1 

k k 

z z 

FIG. 3. Real parts of the dominant eigenvalues of the spec- 
trum as functions of k z : translation (+), rotation (o) and 
meander (•) branches. The kinetics parameters qualitatively 
give small core spirals closer to the 2D meander line in (b) 
than in (a), (a) a = .7, b = .01, and e = .025. The 
value obtained are k\ ~ 2.,uj\ = 1.784 and an = 1.62, 
a± = 0.83 to be compared with a t ~ 1.78i+(-1.63+0.83i) k\ 
and a m ~ -0.157 + 1.94 i + (1.04 + 0.21 i)k% for k z < 1. 
(b) a — .66, b = .01, and e = .025. The value obtained 
are k x ~ 2., ui = 1.756 and ay = 2.44, a± = 0.77 
to be compared with a t — 1.75 i + (—2.41 + 0.75 i)k\ and 
a m ~ -0.071 + 1.89 i + (1.89 + 0.33 i)k'j for k z < 1 . In (b) 
the meander mode is stable at k z — but the meander branch 
becomes unstable at finite k z . 



parameter are those of Fig. ga. When twist is increased 
from zero, the zeroes of Re[a(k z )\ which correspond to 
the translation eigenmodes, move away from k z = 
and stand at k z — ±t w as they should. Most impor- 
tantly, the translation eigenmodes remain local maxima 
of Re[a(k z )]. Beyond a threshold twist (t w ss 0.20) a sec- 
ondary maximum appears for k z close to zero. One of the 
translation branch is shown on Fig. ^a for t w — 0.26 when 
Re[a(k z )] is still negative at this secondary maximum 
(the other translation branch can be deduced from the 
one shown by reflection with respect to the k z = axis). 
When the twist is increased further (past t w w 0.30), 
Re[a(k z )] becomes positive at this secondary maximum. 
The twisted scroll waves then become unstable for a finite 
range of wave numbers as shown on Fig. |]a for t w = 0.4. 
The rotation and meander branches are also shown for 
this value of twist. One sees that the real parts of the two 
meander branches which are superimposed on Fig. ||a for 
t w = have been split by the twist. The maxima around 
k z = ±.48 for t w — have shifted around ±(.48 + t w ) 
but their stability has not been significantly modified. 

Direct numerical simulations show that this twist- 
induced instability leads to a restabilized state in which 
the scroll core itself acquires a helical shape. We can 
thus safely identify the present instability with the 'spro- 
ing' instability of || but again a detailed study of the 
nonlinear state is best deferred to another publication. 

Our results can be partially rationalized on the basis 
of previous works but also show the limitations of some 
previous theoretical suggestions. In the untwisted case, 
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FIG. 4. (a) Real and (b) imaginary parts of the dominant 
eigenmodes for a = .7, b = 0.01, e = 0.025 (same as Fig. 3a) 
and a twist t w = -4: meander(dashed-dotted), rotation (dot- 
ted) and translation (solid) branches. In (a) one of the trans- 
lation branch (bold dashed) is also plotted for a twist t w = .26 
below the instability threshold. 



shown on Fig. 4a remains to be better understood. More 
generally, the precise knowledge of the important linear 
modes and their behavior, should allow the development 
of controlled reduced descriptions of scroll wave nonlin- 
ear dynamics. We also hope that our results and others 
will encourage further experimental studies of scroll wave 
instabilities. 
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it was noted UJ§ that a weak scroll curvature k can be 
transformed away by adding a term —E.Vu on the r.h.s 
of Eq. This relates the curvature-induced motion of a 
weakly curved scroll wave [|l8]j7j to 2D spiral drift in an 
external held E and gives for the small k z behavior of 
the translation branches, 

v±(k z ) = ±iu>i + (— Q!|| ± ia±)k z (5) 

This relation between the translation modes small k z cur- 
vature and the independently measured drift coefficients 
a II and a± |l9[ is very well satisfied by our numerics 
(see the captions of Fig. fj]j|||). In the untwisted case, 
whether translation modes are stabilized or destabilized 
by the introduction of a third dimension is thus precisely 
related to a property of 2D spirals namely, their known 
change of drift direction as kinetics parameters vary poj . 

The analogous question for meander modes is more 
complex than anticipated. A simple normal form pro- 
posed in |8| extends Eq. (||) to the low k z behavior of the 
two complex conjugate meander branches a m ,±(k z ) = 
(J m ,±(0) + (of j| ±ia±)k^. In contrast to Eq. (||), this rela- 
tion is not obeyed by the data of Fig. 0. It even qualita- 
tively disagrees with the data of Fig. Ilia since it predicts 
that the low k z curvature of the translation and mean- 
der branches real parts are of opposite signs contrary to 
what we find. Thus, a simple two-dimensional account of 
the long-wavelength dependence of the meander branch 
remains to be developed. 

For twisted scroll waves, we have found that the mech- 
anism of the 'sproing' instability is a twist-induced de- 
formation of the translation branches. This coupling of 
twist to the translation modes is not present in analyt- 
ical calculations using averaging techniques (l^J^l . We 
will show elsewhere [^lj that it can be analytically cap- 
tured in the weakly excitable limit using the techniques 
of (Tt]]. Nevertheless, the shape of the resulting transla- 
tion branches and specially their double-peak structure 
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